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Abstract 

In this paper an approach is outlined. With this approach some explicit algorithms 
can be applied to solve the initial value problem of n— dimensional damped oscil- 
lators. This approach is based upon following structure: for any non-conservative 
classical mechanical system and arbitrary initial conditions, there exists a con- 
servative system; both systems share one and only one common phase curve; 
and, the value of the Hamiltonian of the conservative system is, up to an additive 
constant, equal to the total energy of the non-conservative system on the afore- 
mentioned phase curve, the constant depending on the initial conditions. A key 
way applying explicit symplectic algorithms to damping oscillators is that by the 
Newton-Laplace principle the nonconservative force can be reasonably assumed 
to be equal to a function of a component of generalized coordinates along a 
phase curve, such that the damping force can be represented as a function analo- 
gous to an elastic restoring force numerically in advance. Two numerical exam- 
ples are given to demonstrate the good characteristics of the algorithms. 

Keywords: Hamiltonian, dissipation, non-conservative system, damping, explicit 
symplectic algorithm 



1. Introduction 

FengjEfl[lll,Marsden||l],Neri0 and Yoshida[7]had developed a series of 
symplectic algorithms for Hamiltonian systems. These algorithms possess some 
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advantages. But it is difficult to apply these algorithms to damping dynamical 
systems, because it has been stated in most classical textbooks that the Hamilto- 
nian formalism focuses on solving conservative problems. Damping phenomena 
is very important in the modeling of dynamical systems, and can not be avoided. 
Our aim is to apply some explicit canonical algorithms to nonlinear damping dy- 
namical systems, which is generated generally by FE-method. These canonical 
algorithms reported in this paper can be readily utilized for computing large-scale 
nonlinear damped dynamical systems. 

Betch[8][9][10] attempted to apply directly some implicit algorithms to damp- 
ing systems. The implicit symplectic algorithms utilized by Betchjsl] possess a 
few good characteristic, e.g. energy-conservation, momentum-consistence, etc... 
In terms of energy-conservation, implicit symplectic algorithms might be better 
than explicit symplectic ones. But explicit symplectic schemes might be more 
suitable for nonlinear problems. 

If one needs to apply symplectic algorithms to a dissipative system, one must 
convert the dissipative system into a Hamiltonian system or find some relationship 
between the dissipative system and a conservative one. 



In the literature [1 1 1 fl . we have stated a proposition describing a relation among 



a damping dynamical system and conservative ones: 

Proposition 1.1. For any non-conservative classical mechanical system and ar- 
bitrary initial condition, there exists a conservative system; both systems sharing 
one and only one common phase curve; and the value of the Hamiltonian of the 
conservative system is equal to the sum of the total energy of the non-conservative 
system on the aforementioned phase curve and a constant depending on the initial 
condition. 

In other words, a dissipative ordinary equation and a conservative equation may 
possess a common particular solution. In the next section, an analytical exam- 
ples are given to explain this proposition. Readers can find the detailed proof of 
Proposition ll.ll in the reference[l 1] 

In the Literature [12] a basic explicit canonical integrator is proposed. Based 
on this basic scheme, Neri^ constructed 4-order explicit canonical integrator, 
and then Yoshida 0] proposed a general method to construct higher order explicit 
symplectic integrator. Utilizing the Proposition ll.il we apply this class of explicit 
canonical integrators to damping dynamical systems. This point will be in detail 
stated in sec. [3] 
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2. One-dimensional Analytical Example 

Consider a special one-dimensional simple mechanical system: 

x + cx = 0, (1) 

where c is a constant. The exact solution of the equation above is 

x=A l +A 2 e- ct : (2) 

where Ai,A2 are constants. Differentiation gives the velocity: 

x=-cA 2 e~ ct . (3) 

From the initial condition xq,xq, we find A\ = xq + xq/c,A 2 = —xq/c. Inverting 
Eq. © yields 

< = -V-^i (4) 
c A 2 

and by substituting into Eq. ©, such we have 

x= -c(x-Ai) (5) 

The dissipative force F in the dissipative system (OQ) is 

F = ex. (6) 

Substituting Eq. © into Eq. ©, the conservative force & is expressed as 

jr = -c 2 (x-Ai); (7) 

Clearly, the conservative force & depends on the initial condition of the dissi- 
pative system (QQ), in other words, an initial condition determines a conservative 
force. Consequently, a new conservative system yields 

x + ^ = 0^x-c 2 (x-A l ) =0. (8) 

The stiffness coefficient in this equation must be negative. One can readily verify 
that the particular solution dU of the dissipative system can satisfy the conserva- 
tive one ([8]). This point agrees with Proposition (11.11) . 
The potential of the conservative system © is 

V = J* [-^(x-Ax)] dx = -yjc 2 + c 2 A lX 



Figure 1 : Relationship between nonconservative system ([T) and conservative one © 



Therefore the Hamiltonian is 

1 c 2 

H = T + V = -p 2 x 2 + c 2 Aix, 

2 2 

where p = x. 

Furthermore, Proposition (ll.lt can be depicted by Fig. |2] The phase flow of 
conservative system © transforms the red area in phase space to the purple area; 
the phase flow of conservative system d8) transforms the red area to the green area. 
The blue curve in Fig. |2] illustrates the common phase curve. If one draws more 
common phase curves and phase flows, the picture will like a flower, the phase 
flow of the nonconservative system likes a pistil and phase flows conservative 
systems like petals. 



3. Modification Symplectic Numerical Schemes 

3.1. Basic Explicit Symplectic Numerical Schemes 

In the paper lfl2h llfjll ifvh a symplectic algorithm based second kind generation 
function was stated 

where the superscript i denotes the i-th time node, q denotes coordinates and p de- 
notes canonical momenta, and H denotes Hamiltonian quantity, H q = dH/dq, H p = 
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dH /dp. If the Hamiltonian is seperable, i.e. H = U(p) +V (q) ,V q = H q , U p = H p , 
then the symplectic scheme© above becomes an explicit symplectic scheme: 



p' +l = p 1 - W^) 

For some nonlinear vibration mechanical system, V q = K(q)q. 
Let us consider an n— dimensional nonlinear oscillator: 

q+Cq+Kq = 0, (11) 

where C denotes a non-linear damping coefficient matrix which depends on q, 
and K denotes a non-linear stiffness matrix which depends on q and consists of 
two parts K = K + K{K is a diagonal matrix). 

In accordance with Proposition 11.11 a conservative mechanical system was 
found associated with the dissipative system (TTTT) in addition to its initial condi- 
tions. Subject to these initial conditions, the dissipative system (fTTI) possesses a 
common phase curve 7 with the conservative system. As in Eq. (|7]), we can con- 
sider that the components of the damping force Cq determine the components of 
a conservative force on the phase curve 7 

cnqi = Pu(qi) ■■■ c\ n q n = p hl (qi) 

\ ••. i (12) 

C n \q\ = P21 (q n ) ■ ■ ■ C nn q n = p nn {q n )- 

For convenience, this conservative force is assumed to be an elastic restoring 
force: 

pll(<?l) = Kn (qi)qi ■■■ pi«(<?i) = K ln (qi)qi 

\ ••• ! (13) 

Pnl (<7l ) = Ki\ {q n )q n ■ ■ ■ Pun (tfn) = Kin (q n )qn- 

In a similar manner, the components of the non-conservative force Kq are 
equal to the components of a conservative force on the phase curve 7 

^1191=^11(^1) •■■ Ki n q n = X\ n {q\) 

: ••• ! (14) 

Kniqi = %2\(q n ) ■ ■ ■ K nn q n = Xnn(qn)- 

The conservative force can likewise be assumed to an elastic restoring force: 

#n(<?i) = ^n(<?iki ••• Xin{qi) = hn {q\)q\ 

\ ••. ; (15) 

Xnl (qi) = k n \ (q n )q n ■ ■ ■ Xnn {qn) = Kin (q n )q n - 
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By an appropriate transformation, an equivalent stiffness matrix K that is diagonal 
in form can be obtained 



i=i 

Consequently, an n-dimensional conservative system is obtained 

q + (K+K)q = 



(16) 



(17) 



which shares the common phase curve y with the n-dimensional damping system 
described by (fTTI) . In this paper, the conservative system is called the 'substitute' 
conservative system. The Lagrangian of Eqs. (fl7T ) is 



(KqYdq- (KqYdq 



with the Hamiltonian 



1 



H=-p T p + 



q 

(Kq) T dq + 



9 _ 



(KqYdq, 



(18) 



(19) 



where is the zero vector, and p = q. Here H in Eq. (fT9~l) is the mechanical energy 
of the conservative system ([T7l ), because /q (Kq) T dq is a potential function such 
that H is independent of the path taken in phase space. 

Subject to a certain initial condition, one need merely to solve the conservative 
system(fT7l). But one must in advance obtain the numerical approximation of the 
matrix K for a time step, such that one can utilize the algorithm (flOl to integrate 
the conservative system (fTTI ) for a time step. One can repeat this process above 
up to the end. In this way one obtains the numerical particular solution of the 
conservative system (fTTI ), which is exactly the numerical particular solution of the 
damping one. The he numerical approximation of the matrix K can be assumed 
as: 



K 



K n 







(20) 



... K nn 

Kj(q/) = Cjiq)/qj l +K jl q l l /qf 

Hence the explicit canonical scheme (ITOT ) can be modified into 

K i j (q/)=c jl qi/q/ + K jl q' l /q/ 
pf+^pf-TlKj+K^q/M) 
q/+ l =q l + tp J '+ 1 



(21) 
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The scheme above is a one order scheme. Furthermore one can construct higher 
order explicit canonical schemes utilizing the method reported in the literatures [Q] . 
Now consider a map from z = z(0) to z' = z( t) : 



z'^([\e ntE e SitF + 0(x n+l ))z, 



(22) 



where 



P 



o 

1 



F 



—(K + K) 




In fact Eq. (|22|) is the succession of the following mappings, 

pJ+l=pJ- Sl W q (qJ) 
q j+1 = q j + riTU p (pj +1 ) ' 



(23) 



In reality the difference between the equations above and Eq.(l2TT) is that the coef- 
ficients S{, r\ before the time step x. In the literature [01 a generalized method to 
determine Si, r ; were given. Therefore, the higher order explicit canonical scheme 
can be represented as: 



K(qJ 







1 

h 
i=l 







K n {q J n ) 
-(K + K) 




KaWa) = I Caiqj/qi+Kaiq'j/q'u 
1=1 



(24) 



4. Numerical Examples 

Two examples will be given to shown this numerical method24l 

4.1. The First Example 

To begin, we consider a Van Der Pol's oscillator 

x + plx(x 2 -\)+x = 0, 



(25) 
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Figure 2: The resonance of the Van Der Pol's oscillator 



where /! = 10. The initial conditions are given by xq = 1, x 
the 4— order explicit symplectic method (|24l with coefficients 



ji=j 4 = [2+(v^+1/v^)]/6, 
n = r 3 = [2 + (^2 + l/^2)]/3, r 2 = 



0. We employee 



5 2 = *3 = [l-(v^+l/v/2)]/6, 
-[2 + ( v / 2 + 1/^)1/3, r 4 = 0, 



and classical explicit 4— order Runge-Kutta method to compute the resonance of 
the Van Der Pol's oscillator (|26l respectively, then employ a same method to in- 
tegrate the results to the total energy, which is the sum of the mechanical energy 
and the work done by damping forces in the system (|25l ). The both methods are 
run with the same step size X = 0.01. The resonance is shown in Fig. |2l and the 
total energy is shown in Fig. |3] 

It is aparent from Fig. [3] that the explicit symplectic method (|24l) has qual- 
itatively different behavior to the Runge-Kutta method. The energy divergence 
between the explicit symplectic method and the exact solution is smaller than 
that between Runge-Kutta method and the exact solution. The energy divergence 
between the explicit symplectic method and Runge-Kutta method increases with 
the time evolution. Due to the increasement of the energy, the phase difference 
between both the results in Fig. |2]increases also with the time evolution. 
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Explicit Symplectic 
Runge-Kutta 
Exact solution 



20 30 40 50 60 70 80 90 

t 



Figure 3: The total energy of the Van Der Pol's oscillator 



4.2. The Second Example 

In the second example, we consider a 2— dimensional damped nonlinear Duff- 
ing oscillator 

2q\ +0.\qi + (2 + 0.\qj)qi +42 = 

3q 2 + 0.2q 2 + qi + (2 + 0.2^)^ = 0, 1 ' 

with the initial conditions q\ =0, q 2 = 0, q\ = 0, q 2 = 1. The program of the 
both methods with the step size % = 0.01 are carried out to simulate Eq. (l26l ). The 
resonance is shown in Fig. |2] the numerical solution of the total energy is shown 
in FigjU 

There is only tiny difference between resonance results of the two methods, 
correspondingly, the difference among the total energy obtained by the numerical 
methods and anlytical methods is very tiny. As numerical examples in the other 



literatures II 1311 . that explicit Runge-Kutta method must cause numerical pseudo 



dissipation which might be positive or negative. The difference between our nu- 



merical examples and the examples in the literature II 13] 
examples and the mechanical energy in their examples^ 



is the total energy in our 



Fig. 6.1 in the literature[13] 
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Explicit Symplectic 
Runge-Kutta 




20 40 60 80 100 120 140 160 180 200 220 

t 

Figure 4: The 1-th displacement of the damped Duffing oscillator 




Explicit Symplectic — h- 
Runge-Kutta 
Exact solution - - - -X - 



20 40 60 80 100 120 140 160 180 200 

t 



Figure 5: Total energy of the damped dissipative oscillator 
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5. Conclusions 



We have introduced a class of explicit symplectic algorithms to dissipative me- 
chanical systems successfully, by changing these algorithms into the scheme .(l24l). 
Because the algorithms (|24l) are explicit and possess good energy preserving char- 
acteristics, the explicit symplectic algorithms (1241 is quite suitable for long term 
integration of arbitrary dimensional nonlinear dissipative mechanical systems. 
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